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The mid-Pliocene (-3 to 3.3 Ma ago), is a period of sustained global warmth in comparison to the late 
Quaternary (0 to -1 Ma ago), and has potential to inform predictions of long-term future climate change. 
However, given that several processes potentially contributed, relatively little is understood about the rea- 
sons for the observed warmth, or the associated polar amplification. Here, using a modelling approach and 
a novel factorisation method, we assess the relative contributions to mid-Pliocene warmth from: elevated 
CO 2 , lowered orography, and vegetation and ice sheet changes. The results show that on a global scale, the 
largest contributor to mid-Pliocene warmth is elevated C0 2 . However, in terms of polar amplification, 
changes to ice sheets contribute significantly in the Southern Hemisphere, and orographic changes contribute 
significantly in the Northern Hemisphere. We also carry out an energy balance analysis which indicates that 
that on a global scale, surface albedo and atmospheric emmissivity changes dominate over cloud changes. We 
investigate the sensitivity of our results to uncertainties in the prescribed C0 2 and orographic changes, to 
derive uncertainty ranges for the various contributing processes. 

© 2012 Elsevier B.V. All rights reserved. 


1. Introduction 

The most recent palaeoclimate reconstructions (Dowsett et al., 
2009) suggest that during warm ‘interglacials' of the Pliocene epoch 
(-5.3 to 2.6 Ma), global annual mean sea surface temperatures 
were 2 to 3 °C higher than the pre-industrial era. During these 
warm interglacials sea levels were higher than today (estimated to be 
10 to 30 + m) meaning that global ice volume was reduced (e.g. 
Dowsett and Cronin, 1990; Dwyer and Chandler, 2009; Naish and 
Wilson, 2009). There were large fluctuations in ice cover on Greenland 
and West Antarctica, and during the interglacials they were probably 
largely free of ice (Dolan et al., 2011; Hill et al., 2010; Lunt et al., 
2008; Pollard and DeConto, 2009). Some ice may also have been lost 
from around the margins of East Antarctica especially in the Aurora 
and Wilkes sub-glacial basins (Hill et al., 2007). Coniferous forests 
replaced tundra in the high latitudes of the Northern Hemisphere 
(Salzmann et al., 2008), and the Arctic Ocean may have been seasonally 
free of sea-ice (e.g. Cronin et al., 1993). The most recent estimates of 
Pliocene atmospheric C0 2 concentrations range between 280 and 
450 ppmv (Pagani et al„ 2010; Seki et al., 2010). The Mid-Piacenzian 
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Warm Period (henceforth ‘mid-Pliocene’; 3.26 to 3.025 Ma BP; time- 
scale of Lisiecki and Raymo (2005)) is a particularly well documented 
interval of warmth during the Pliocene, with global data sets of 
multi-proxy sea surface temperatures, bottom water temperatures, 
vegetation cover, topography and ice volume readily available as 
boundary conditions and/or evaluation datasets for global climate 
models (Dowsett et al., 2010b; Haywood et al„ 2010). 

Many parallels have been drawn between the apparent similari- 
ties in climate between warm intervals of the Pliocene and the end 
of the 21st Century, particularly in terms of (relative to pre- 
industrial) (a) the change in annual mean global temperature 
(Haywood et al., 2000a; Jansen et al., 2007), (b) changing meridional 
surface temperature profiles showing a strong polar amplification of 
the warming (Dowsett et al., 1992; Robinson, 2009), (c) changing 
precipitation patterns and storm tracks (Haywood et al., 2000b) and 
even (d) Hurricane intensity and ENSO-event frequency/extra tropi- 
cal teleconnections (Bonham et al., 2009; Fedorov et al., 2010; 
Scroxton et al„ 201 1 ; Watanabe et al„ 201 1 ). This attraction is made 
more intense by the fact the continents had essentially reached 
their modern position, and due to its relative youth, geologically 
speaking, inferences about the environmental tolerances of many of 
the biological proxies used to reconstruct Pliocene environments 
and climates can be made with far greater confidence than further 
back in Earth history (Dowsett and Poore, 1996; Salzmann et al., 
2008). 
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As such, it is particularly important to understand why the mid- 
Pliocene was warmer than pre-industrial. Up until now, the most 
comprehensive attempt to answer this question was carried out by 
Haywood and Valdes (2004), henceforth H&V04. Using the UK Met 
Office coupled atmosphere-ocean General Circulation model, 
HadCM3, they carried out a model simulation of the mid-Pliocene, 
and compared it to a pre-industrial simulation. They found a global 
mean surface air temperature difference of 3.1 °C. From the assumed 
C0 2 radiative forcing in the model and consideration of top-of-the- 
atmosphere radiative fluxes, they partitioned the causes of this tem- 
perature difference between C0 2 (l,9Wm‘ 2 ), surface albedo 
(2.3 Wm~ 2 ) and cloud cover (1.8 Wm -2 ) changes. They further 
partitioned the surface albedo component between land ice and 
snow (55%) and sea ice (45%) changes. From interrogating the 
ocean stream function and net heat transports, they also concluded 
that ocean circulation changes did not lead to significant surface tem- 
perature warming. Given the considerable computational constraints 
at the time (the 300 yr simulation took 9 months to complete), the 
H&V04 study contributed significantly to our understanding of the 
causes of mid-Pliocene warmth. However, the fact that further sensi- 
tivity studies could not be carried out meant that cause and effect was 
not easily partitioned. For example, the albedo change due to sea ice 
was itself a result of the imposed C0 2 (and orography, and vegetation, 
and land-ice) changes. Similarly for clouds — some of the cloud 
changes would be due to the land ice (and other) changes. In this 
paper we address this issue, by describing a new methodology for a 
robust, self-consistent partitioning of climate change between several 
causal factors. We then apply it to the warm periods of the mid- 
Pliocene, resulting in a partitioning of temperature changes between 
changes in the prescribed C0 2 , orography, vegetation and ice sheet 
boundary conditions. We also carry out an analysis of the pre- 
industrial and mid-Pliocene results using an energy balance method 
described by Heinemann et al. (2009). 


2. Experimental design 

2.1. Model description — HadCM3 

All the General Circulation Model (GCM) simulations described in 
this paper are carried out using the UK Met Office coupled ocean- 
atmosphere GCM HadCM3, version 4.5 (Gordon et al., 2000). The res- 
olution of the atmospheric and land components is 3.75° in longitude 
by 2.5° in latitude, with 19 vertical levels in the atmosphere. The 
resolution of the ocean model is 1.25° by 1.25° with 20 levels in the 
vertical. Parameterisations include the radiation scheme of Edwards 
and Slingo (1996), the convection scheme of Gregory et al. (1997), 
and the MOSES-1 land-surface scheme, whose representation of evap- 
oration includes the dependence of stomatal resistance on tempera- 
ture, vapour pressure and C0 2 concentration (Cox et al., 1999). The 
ocean model uses the Gent and McWilliams (1990) mixing scheme. 
There is no explicit horizontal tracer diffusion in the model. The hori- 
zontal resolution allows the use of a smaller coefficient of horizontal 
momentum viscosity leading to an improved simulation of ocean ve- 
locities compared to earlier versions of the model. The sea ice model 
uses a simple thermodynamic scheme and contains parameterisations 
of ice concentration (Hibler, 1979) and ice drift and leads (Cattle and 
Crossley, 1995). In simulations of the present-day climate, the model 
has been shown to simulate SST in good agreement with modern ob- 
servations, without the need for flux corrections (Gregory and 
Mitchell, 1997). Future climate predictions from the model were pre- 
sented in the latest 1PCC report (Solomon et al., 2007), and it has been 
used in the Palaeoclimate Modelling Intercomparison Project to simu- 
late Last Glacial Maximum and Mid-Holocene climates (Braconnot et 
al., 2007). The model will also be used in the forthcoming PlioMIP 
project (Haywood et al., 2010, 2011b). 


2.2. Boundary conditions 

The PRISM project (http://geology.er.usgs.gov/eespteam/prism/) 
has as its main aim the characterisation of the palaeoenvironment 
of the mid-Pliocene warm period (3.26-3.025 Ma) on a global scale. 
In this paper, we simulate the mid-Pliocene climate by making use 
of the PRISM2 reconstruction of orography, vegetation, and ice sheet 
extent (Dowsett, 2007; Dowsett et al., 1999), which are described 
below. 

The PRISM2 orography reconstruction was based on palaeobotani- 
cal evidence suggesting that the East African rift areas were 500 m 
higher during the mid-Pliocene relative to today (Thompson and 
Fleming, 1996). In contrast, palaeoelevation of the Western Cordillera 
of North America and northern South America was reduced by 50%. 
Large elevation differences are noted in both Greenland and Antarcti- 
ca due to significant removal of continental ice (Dowsett, 2007; 
Dowsett et al., 1994). PRISM2 land ice distribution and volume was 
closely associated with sea level estimates from several sources (see 
Dowsett, 2007), which indicate a eustatic sea level rise of around 
25 m compared to modern. These estimates have recently been con- 
firmed by independent studies based on the depth palaeoecology of 
foram assemblages from New Zealand (Naish and Wilson, 2009) 
and benthic Mg/Ca and oxygen isotopes (Dwyer and Chandler, 
2009). Antarctic ice distribution was based upon a modelled stable 
ice sheet configuration (see Dowsett et al., 1999), strongly con- 
strained by the sea-level reconstructions. The PRISM2 vegetation re- 
construction (Dowsett et al., 1999) was compiled from fossil pollen 
and plant macrofossil data from 74 sites covering all continents. 
PRISM2 vegetation is identical to PRISM1 (see Thompson and 
Fleming, 1996). PRISM2 uses seven land cover categories (desert, tun- 
dra, grassland, deciduous forest, coniferous forest, rainforest, and land 
ice) that are a simplification of the 22 land cover types of Matthews 
(1985). From the PRISM2 vegetation, orography, and ice-sheet extent, 
we derive all the boundary conditions necessary to run the GCM in 
mid-Pliocene mode (a total of 23 variables different to those of the 
pre-industrial, such as heat capacity of the soil, albedo, moisture hold- 
ing capacity etc.). 

Since the development of the PRISM2 dataset, the USGS have now 
released an updated version — PRISM3 (Dowsett et al., 2010a, 2010b). 
We use the PR1SM2 dataset; firstly, to maintain consistency with pre- 
vious modelling studies, in particular H&V04 and Lunt et al. (2010a); 
secondly, the mid-Pliocene simulation with PR1SM2 boundary condi- 
tions has been spun up for a total of over 1 000 yrs, which is consider- 
ably more than could be achieved with new boundary conditions in a 
reasonable timeframe. In Section 4.1 we discuss the implications for 
this study of using PRISM2 compared to PRISM3 boundary conditions. 

2.3. Factorisation methodology 

The primary aim of this study is to assess the relative importance 
of various boundary condition changes which contribute to mid- 
Pliocene warmth. Therefore, we are aiming to partition the total 
mid-Pliocene warming, AT, into four components, each due to the 
change in one of the boundary conditions C0 2 , orography, ice sheet, 
and vegetation. The assumption here is that other palaeogeographic 
changes not currently captured by the PRISM dataset, such as soils 
or lakes, have a negligible impact on the global mean temperature 
change. 

AT = dT C02 + dT orog + dT ice + dT veg (1) 

‘Factor separation’ techniques (e.g. Stein and Alpert, 1993) can be 
used to determine these components of the mid-Pliocene surface air 
temperature change dT C o 2 , dT orog , dT veg , and dT ice . Typically, this in- 
volves carrying out an ensemble of GCM simulations with various 


130 


D.J. Lunt et al. / Earth and Planetary Science Letters 321-322 (2012) 128-138 


combinations of boundary conditions. Here we present a new factor- 
isation methodology, which we believe improves on previous work. 

We name a GCM simulation which has boundary conditions x and 
y modified from pre-industrial to mid-Pliocene as E xy . The four 
boundary conditions considered are atmospheric C0 2 (c), orography 
(o), vegetation (v), and ice sheets (i). Thus, a pre-industrial simula- 
tion is E, a mid-Pliocene simulation is E ociv , and e.g. a simulation 
with pre-industrial ice sheets and vegetation but mid-Pliocene orog- 
raphy and C0 2 is E oc . The corresponding surface air temperature dis- 
tributions in these simulations we name T, T ociv , and T oc respectively. 

For simplicity, we first describe our factorisation methodology by 
considering a simpler example, where only two boundary conditions 
(C0 2 and orography) are changed instead of four. The simplest factor 
separation technique is the incremental application of the boundary 
conditions. For our simplified example, this could involve an ensem- 
ble of 3 GCM simulations: £, £ c , and £ oc . The total temperature anom- 
aly, AT (equal to T oc —T in this simplified example), could be 
separated into 2 components: 

dT C02 =T c -T 

dT omg = T oc -T c , W 


This method, illustrated in Fig. la, has been used extensively in the 
climate literature (e.g., for the LGM see Broccoli and Manabe, 1987; 
von Deimling et al., 2006). It has the advantage that a limited number 
of simulations (N + 1, where N is the number of processes investigat- 
ed) need be carried out. It has the disadvantage that it results in a 
non-unique solution: one could equally define 

dTco 2 =T 0C —T 0 

dToro g = T 0 -T, ( j 

which, due to non-linearities would in general result in a different 
partitioning. 

Stein and Alpert (1993) (henceforth S&A93) recognised this and 
instead suggested that, considering the temperature response as a 
continuous function of two variables (in our simplified example orog- 
raphy and C0 2 ), and carrying out a Taylor expansion about the con- 
trol climate, one can write 

AT = ACO - 1 + . ^ Aorog + nonlinear terms. (4) 

oC0 2 aorog 


They suggested that the nonlinear terms could be considered as 
‘synergy’, S, between the two forcing variables, and that the partial 
derivatives be estimated from the GCM simulations relative to the 
control, so that 


This method, illustrated in Fig. lb, has been used in several previ- 
ous studies (e.g. for the mid-Holocene and LGM see Jahn et al., 2005; 
Wohlfahrt et al., 2004). It has the advantage that it takes into account 
the non-linear interactions between the different boundary condi- 
tions. However, it requires a larger number of simulations (2 N ) than 
the linear approach. Perhaps more importantly, it has the problem 
that it is not symmetric: one could equally carry out the Taylor expan- 
sion about the perturbed climate, and write 

~dT C o 2 = T 0 —T 0C 

-dT orog = T C ~T 0C (6) 

-S = T-T 0 -T C + T 0C 

i.e. it would in general give a different answer if one asked “why is the 
mid-Pliocene warmer than pre-industrial" than if one asked “why is 
the pre-industrial cooler than the mid-Pliocene” (although the syner- 
gy term, S, would have the same magnitude in both cases). 

In order to obtain a symmetric and unique factorisation, we in- 
stead estimate the partial derivatives in Eq. (4) with their average 
values over the domain considered, and write for our simplified case: 

dT C o 2 =\((T c -T) + (T oc -T 0 )) 
dTomg = r;((T 0 -T) + (T oc -T c )). 

This is equivalent to averaging the two different formulations of the 
S&A93 approach in Eqs. (5) and (6). An alternative, but identical, inter- 
pretation is that our technique uses the S&A93 formulation of Eq. (5) 
but attributes the synergy term, S, equally between the two forcings: 

dTco 2 = T c —T + S/2 

dT 0mg = T 0 -T + S/2 (8) 

(S = T 0C -T 0 -T C + T). 

It is also equivalent to averaging the two linear formulations in 
Eqs. (2) and (3). 

Our formulation has the advantage that it takes into account non- 
linear interactions, and is symmetric. In common with the S8 jA 93 
approach, it requires 2 N GCM simulations, and so is more computa- 
tionally demanding than the linear approach. 

For our mid-Pliocene study, where we actually have 4 variables 
(C0 2 , orography, vegetation, and ice sheets), this would require 
2 4 = 16 simulations. The factorisation would be as follows: 

dT C o 2 «g((T e -T) + (T oc —T 0 ) + (T ic — T,) + (T vc -T v ) (9) 

+(T ocv —T ov ) + (T oci — T oi ) + (Jciv—Ti V ) + (T odv —T oiv )), 


dT C02 =T c -T 
dT 0 m g — T 0 —T 
S=T 0C -T 0 -T C + T. 


( 5 ) 


dT orog - g ((T 0 —T) + (T co —T c ) + (T jo —Tj) + (T vo — T v ) ^q-j 

+{T C0V —T CV ) + (T coi —T d ) + (T oiv —T iv ) + ( T c0 j v — T clv )), 


(a) 


Toc 


Tc 


dTco2=Tc-T 


C02 


(b) 


dTorog=Toc-Tc 


dTorog=To-Tc 


Tc 


dTco2=Tc-T 


C02 


(c) 


S=Toc-To-Tc+T 


dTorog=l/2((To-T)+( Toc-Tc )) 
dTco2 = l/2((Tc- T)+(Toc—To)) 


C02 


Fig. 1 . Factor separation for a function of two variables— in this case C0 2 and orography, (a) is the linear approach (Eq. 2), (b) is the Stein and Alpert (1993) approach (Eq. 5), and 
(c) is our approach (Eq. 7 or Eq. 8). 
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Tocv 


Tocvi 


dTveg=l/2((Tocv-Toc) 

dTice— 1/2( (Toci—Toc 


+(Tocv 
)ff (Tocvi 


-Toci)) 

Tocv)) 


the modern ice. In other words, it is not well defined whether the 
albedo-induced warming associated with reduced ice sheets during 
the Pliocene is due to the reduction of ice per se, or due to the vegeta- 
tion which replaces it. Here, we make the decision to attribute this 
warming to the vegetation that replaces it. As such, both simulations 
E oci and £ ocv have the albedo of ice in regions which are ice-free in the 
Pliocene but have ice in the modern (Fig. 3). 


To 


Toe 


dTorog=l/2((To-T)+( 

dTco2=l/2((Tc-T)+(T( 


Toc-Tc)) 
oc -To)) 


C02 


Fig. 2. Factor separation used in our study for two functions of two variables each — in 
this case C0 2 , orography, vegetation, and ice (Eq. 13). 


dT veg — g ((T v — T) + (T„~r c ) + (T iv — Ti) + ( T 0V —T 0 ) 

+(^0,0“ ^co) + devi Td) + (T vi0 -TJ + (Tn™,— T do )), 


dTj ce 


l«Ti-T) + (T ci -T c ) + (T vi -T v ) + (T 0i -T 0 ) 
+( T cio~ T co) + (Jci V — T cv ) + (T ivo -T vo ) + (T civo — !„„)). 


( 12 ) 


2.4. Mid-Pliocene model-data comparison 

Before presenting and discussing our results, it is first important to 
have some confidence that the mid-Pliocene simulation, E ociv , is con- 
sistent with observations of that period. 

The SSTs in our mid-Pliocene simulation were evaluated relative 
to reconstructions of mid-Pliocene SST in Lunt et al. (2010a). They 
showed that the global mean SST change, mid-Pliocene minus pre- 
industrial, was well simulated (1.83 °C in the model and 1.67 °C in 
the observations). However, they also found that the latitudinal 
distribution of temperature change was not well simulated (their 
Fig. 3c); the modelled mid-Pliocene warming being too great in the 
tropics and too small towards the poles. These discrepancies were 
investigated and discussed further in Dowsett et al. (2011). 

A model-data comparison for the terrestrial climate, using a data- 
base of Pliocene palaeobotanical data (Salzmann et al., 2008, 2009) 
was presented in the Supplementary Information of Lunt et al. 
(2010a) They found a fair agreement between E ociv and the data on 
a global scale, with significantly improved skill at high latitudes in 
the Eodv simulation compared with the pre-industrial E simulation. 

3. Results 


Given the computational expense of carrying out 16 fully-coupled 
GCM simulations, we choose instead to consider C0 2 /orography, and 
vegetation/ice sheets separately, and carry out two N = 2 factor sepa- 
rations (as in Eq. 13), requiring only 7 simulations (illustrated in 
Fig. 2). 

dT co 2 =^((T c -T) + (T 0C -T 0 )), 

dT omg = h(T„-T) + (T oc -T c )), 

z (13) 

dT veg = 2 ((T 0 cv — T 0 c) + (T OC iv~T oci )), 
dTice = 2 ((Toci — Toc) + bociv ^ocv))- 

This factorisation is more computationally efficient than the full 
factorisation in Eqs. (9-12), but is not fully symmetric. 

Five of these simulations ( E , E 0 , E c , E oc , £ ociv ) were used in the study 
of Lunt et al. (2010a) in the context of deriving estimates of Earth sys- 
tem sensitivity, and the orography and snow-free surface albedo of 
these simulations are shown in their Table 1 of their Supplementary 
Information. The orography and snow-free albedo (an indicator of 
the land ice and vegetation distributions) for the 2 new simulations 
( E oci , E 0 cv), along with those for E and £ ociv for comparison, are 
shown in Fig. 3. It is worth noting that because the ice sheets and veg- 
etation are mutually exclusive in any one model grid cell, it is not pos- 
sible to uniquely define boundary conditions for simulations £ oci and 
E ocv . For the simulation with modern vegetation but Pliocene ice 
sheets (£ 0 d). in the regions which are ice sheet-free in the Pliocene 
but have ice sheets in the modern (e.g. the West Antarctic peninsula), 
it is not clear what albedo should be prescribed as there is no modern 
vegetation defined in these regions. Similarly, for the simulation with 
modern ice but Pliocene vegetation (£ ocv ), in the same regions it is 
unclear whether to use the albedo of the Pliocene vegetation or of 


The temperature changes due to the C0 2 (dT C c> 2 ), orography 
(dTo m g), vegetation ( dT veg ) and ice sheet (dT ice ) boundary condition 
changes, as calculated from Eq. (13), as well as the total change, AT, 
are illustrated in Fig. 4. As a global average, of the total mid- 
Pliocene 3.3 °C temperature change, 1.6 °C (48%) is from the C0 2 
(dT C o 2 ), 0.7 °C (21%) is from the orography (dT orog ), 0.7 °C (21%) is 
from the vegetation ( dT veg ), and 0.3 °C (10%) is from the ice sheets 

(dTice). 

dT C o 2 (Fig. 4b) represents the temperature change due to C0 2 
alone. It shares much in common with similar (C0 2 doubling as op- 
posed to 280-400 ppmv here) results presented in the most recent 
report of the IPCC (Solomon et al., 2007). For example, there is polar 
amplification due to snow and sea ice feedbacks, and greater temper- 
ature change on land compared to ocean due to reduced latent cool- 
ing and lower heat capacity. The North Atlantic shows reduced 
temperature increase due to ocean mixing and reduced northward 
heat transport in the Atlantic due to an increase in the intensity of 
the hydrological cycle. The increase of 1.6 °C implies a climate sensi- 
tivity due to a doubling of C0 2 of -3.2 °C, which is close to the middle 
of the IPCC range (Solomon et al„ 2007). dT orog (Fig. 4c) highlights the 
local lapse-rate warming effect of the lower mid-Pliocene Rocky 
Mountain range. There is also a cooling to the west of the mid- 
Pliocene Canadian Rockies, associated with reduced precipitation 
and cloud cover, due to reduced ascent over the mountain range. 
There is a significant non-local effect of the lower Rockies — there is 
a large Arctic warming, in particular in the Barents Sea, which is am- 
plified by reduced sea ice cover. This is due to a modification of the 
Rossby wave pattern, which is more zonally symmetric with the 
lower Rockies, indicated by a reduced trough over Greenland in the 
500 mbar geopotential height field, consistent with previous work 
(e.g. Foster et al., 2010; Kutzbach et al., 1989). Veiy localised cooling 
associated with topographic effects are seen in the Andes, Himalayas, 
and East African rift valley regions. The surface ocean warming east of 
Japan is consistent with previous work showing this to be a region 
sensitive to orographic change in this model (Lunt et al., 2010b). 
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Fig. 3. Orography and snow-free albedo for the E, E oci , E ocv , and E ociv GCM simulations. For equivalent figures of the other GCM simulations ( E 0 , E c , and E oc ), see Table 1 of Supple- 
mentary Information of Lunt et al. (2010a). 
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Fig. 4. (a) Simulated annual mean surface air temperature change, mid-Pliocene minus pre-industrial, AT. (b-e) Surface air temperature changes due to (b) C0 2 ( dT C o 2 ), (c) orog- 
raphy (dT„ro S ), (d) vegetation (dT veg ), and (e) ice (dT,„); as calculated from Eq. (13). 


dT ve g (Fig. 4c) shows that the largest vegetation-related temperature 
changes are in the Canadian Arctic, in particular Greenland (change 
from ice sheet to boreal forest), the Canadian archipelago (change 
from bare soil and glaciers to boreal forest), and Siberia (change 
from bare soil to boreal forest). This warming can be attributed to 
the relatively low albedo of boreal forest in the model, even when 
there is snow-cover on the ground. There are also large changes in 
the tropics, in particular in the Arabian peninsula, where the 


PRISM2 reconstruction indicates a shift from desert to grassland veg- 
etation (based on pollen data (Van Campo, 1991)), resulting in a 
lower albedo in the mid-Pliocene than in the modern (see Fig. 3). 
Some of the temperature changes attributed to vegetation will also 
be due to modifications to the roughness length, potential evapo- 
transpiration, and other vegetation-specific model parameters. dT ice 
(Fig. 4d) shows warming in Greenland and parts of Antarctica due 
to a combination of lapse-rate, due to a lower mid-Pliocene ice 
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sheet height, and albedo, due to the less reflective mid-Pliocene sur- 
face. The regions of Antarctic cooling are due to the fact that the 
PRISM ice sheet is higher in the Pliocene than in the modern in 
these regions. This is consistent with increased precipitation in the in- 
terior of the East Antarctic ice sheet in the warmer climate, and with 
modelled predictions for the future evolution of the Antarctic ice 
sheet under greenhouse gas forcing (e.g. Huybrechts et al., 2004). 
The cooling in the Barents Sea is also consistent with previous work 
investigating the climatic effects of the removal of the Greenland ice 
sheet (Lunt et al., 2004; Toniazzo et al., 2004). However, apart from 
in this region, the signal due to the removal of the ice is very localised. 

The results also allow us to ascertain the contribution to polar am- 
plification of the four factors. We define polar amplification in this 
case to be any warming in the polar regions which is greater than 
the global mean warming. Fig. 5(a) shows the same results as in 
Fig. 4, but as zonal means. It is clear that the polar amplification in 
the Southern Hemisphere is due primarily to the ice sheet changes, 
whereas in the Northern Hemisphere it is due primarily to a combina- 
tion of C0 2 and orography changes, with some contribution from veg- 
etation around 60-70°N. Figs. 1-3 in Supplementary Information 
illustrate the seasonality of the factorisation and polar amplification. 
It is clear that in the Northern Hemisphere, the polar amplification 
is dominated by an autumn and winter signal; in JJA there is 
almost no Northern Hemisphere polar amplification. In the Southern 




Fig. 5. (a) Zonal annual mean surface air temperature changes due to C0 2 (dT C o 2 ), 
orography (dT orog ), vegetation (dT veg ), and ice ( dT ice ) [°C], (b) Zonal annual mean surface tem- 
perature changes due to various components as derived from the energy balance analysis de- 
scribed in Heinemann et al. (2009). Also shown is the total change from the analysis (EBM) 
and from the climate model itself (GCM). 


Hemisphere the seasonality is much more muted. These features are 
consistent with sea-ice and snow being the main causes of the 
seasonality. 

It is interesting to assess the linearity of the climate system to 
these changes in boundary conditions. For example, to what extent 
does the temperature response of the system to a C0 2 change depend 
on the climate base state. Or, in other terms, how large is the ‘synergy' 
term (S in Eq. 5) in the S&A93 formulation? Fig. 6 shows the two 
terms (T c — T and T oc — T„ ) which make up dT C o 2 in Eq. (8), and the 
difference between them (S). The non-linearity is small compared 
to the temperature change itself, showing that in this case, the tem- 
perature response to an increase in C0 2 is largely independent of 
the orographic configuration. Similarly, the vegetation and ice sheet 
changes exhibit relatively small non-linearity (not shown). This im- 
plies that in this case, similar results could be obtained with a simple 
linear factorisation. However, it is not possible to know this a priori. 
The subtle non-linearities of the response of the system to changes 
in C0 2 alone are discussed in more detail in Haywood et al. (2011a). 

It is also instructive to compare our results with those of H&V04. 
Our mid-Pliocene simulation differs from that of H&V04 for two rea- 
sons. Firstly, our simulation is a continuation of that of H&V04, and so 
is further spun-up and closer to equilibrium. Secondly, our simulation 
has been carried out over a number of ‘real-world’ years, and over this 
time has been migrated across several computers and Fortran com- 
pilers. Both hardware and compiler changes can affect the mean equi- 
librium climate of a model, due at least in part to non-standard 
programming practice, for example multiple ‘data’ statements in For- 
tran subroutines (Steenman-Clark, 2009). Fig. 7a shows the differ- 
ence in mid-Pliocene surface air temperature between our 
simulation and that of H&V04, and Fig. 7b shows the difference in 
mid-Pliocene surface air temperature anomaly, mid-Pliocene minus 
pre-industrial, between our simulation and that of H&V04. Our mid- 
Pliocene simulation is significantly cooler than that of H&V04 
( — 0.8 °C in the global annual mean), but the difference in anomalies 
is smaller (0.3 °C). Examination of the temporal evolution of these 
differences indicates that the effect of hardware and compiler change 
is more important than the effect of increased spinup time. This un- 
derlines the importance of always carrying out sensitivity simulations 
on the same machine, and with the same compiler, as any control 
simulation. 

As stated in the Introduction, H&V04 estimated the contributions to 
mid-Pliocene warmth by considering aspects of the global energy bal- 
ance. Heinemann et al. (2009), in the context of the Eocene, present a 
different method of energy-balance analysis which includes a meridio- 
nal analysis. Here, we use the method of Heinemann et al. (2009) to an- 
alyse our mid-Pliocene (E ociv ) and pre-industrial (£) simulations. The 
method gives latitudinal distributions of the contribution to the surface 
temperature change, E ociv — E, of: (a) emissivity changes due to changes 
in greenhouse gases, (b) emissivity changes due to changes in clouds, 
(c) albedo changes due to changes in the planetary surface, (d) albedo 
changes due to changes in clouds, and (e) heat transport changes. This 
latitudinal partitioning is shown in Fig. 5(b). The first thing to note is 
that this approach is based on zonal and seasonal means, and as such 
the total surface temperature change is slightly underestimated by the 
energy balance approach (compare the green line with the black line 
in Fig. 5(b)). On a global scale, the contribution of heat transports to 
the total change is by definition zero, but in the Northern Hemisphere 
there is a small positive contribution at high latitudes and a small nega- 
tive contribution at low latitudes, consistent with a slight increase in 
poleward heat transport. The global mean contribution of clouds (both 
albedo and emmissivity effects) is relatively small, but in the short- 
wave this results from a cancellation of a positive contribution in the 
tropics and a negative contribution at mid-high latitudes. Changes in 
emmissivity (due to the increase in greenhouse gas from 280 to 
400 ppmv, and the associated water vapour forcing) contribute 61% of 
the total surface temperature change, with greatest contribution in 
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mid-high latitudes. Surface albedo changes contribute 44%, due almost 
entirely to mid-high latitude changes; in the tropics the change in albe- 
do contributes very little. Overall it can be seen that surface albedo and 
direct greenhouse-gas forcing are the greatest contributors to the total 
change, with the greenhouse gas forcing dominating in low latitudes, 
and the surface albedo changes dominating at mid-high latitudes. The 
polar amplification is significantly dampened by changes in short- 
wave cloud forcing. It should be noted that cloud processes are amongst 
the most uncertain in GCMs, and so these results are likely to be model 
dependent. 

4. Discussion 

Here we discuss some of the assumptions in this work, including 
quantitative estimates of how some of these assumptions could affect 
our results. 

4.1. Palaeoenvironmental boundary conditions 

In Section 2.2 we describe why we use the PRISM2 boundary con- 
ditions as opposed to the PR1SM3 boundary conditions. The most 
significant effect of this is likely related to the different orography 
dataset in PRISM3 compared to PRISM2 (the ice sheets, although dif- 
ferent, are similar in extent and height, and the PRISM3 vegetation is 
based on an extended dataset which includes PRISM2 as a subset). 
PRISM3 orography is based on the reconstruction of Markwick 
(2007). It differs from PR1SM2 mainly in the high Eurasian latitudes 
and the Himalayas where the geological evidence is inconclusive 
and debated (e.g. Rowley and Garzione, 2007; Spicer et al., 2003). 
The Markwick (2007) reconstruction is actually much closer to mod- 
ern than that of PRISM2. Therefore, using modern orography instead 
of PR1SM2 provides an end-member approximation for the uncertain- 
ty in our results. In this case, given the linearity of the system 
highlighted in Section 4, the total mid-Pliocene temperature change 
can be approximated by: 

Ar „oorog = AT _ dTomg = ^ + ^ + dTke ( 14) 

which is 2.6 °C. Then, the partitioning (Table 1) is 1.6 °C (61%) from 
the C0 2 (dTcoJ, 0.7 °C (27%) is from the vegetation (dT veg ), and 
0.3 °C (13%) from the ice sheets (dT ice ). 

There is no information given in either PRISM2 or PRISM3 on pos- 
sible bathymetric differences between the mid-Pliocene and present. 
As such, we use modern bathymetry in the simulations presented 
here. However, geophysical records of mantle temperature beneath 
the North Atlantic indicate that the Greenland-Scotland ridge was 
about 300 m lower in the Pliocene than modern (Robinson et al., 
2011). A recent modelling study (Robinson et al., 2011) has shown 
that, although this has negligible effect on the global mean tempera- 
ture, it could lead to increased polar warmth (greater than 5 °C) in the 


mid-Pliocene due to increased oceanic northward heat transport in 
the North Atlantic. This has the effect of bringing the modelled SSTs 
in the mid-Pliocene E odv simulation into better agreement with the 
PRISM3 proxy estimates in this region. 


4.2. Mid-Pliocene C0 2 

Mid-Pliocene atmospheric C0 2 has been reconstructed by a varie- 
ty of proxies. A value of 400 ppmv has been used in this and several 
other previous modelling studies of the mid-Pliocene climate (includ- 
ing H&V04), but there are uncertainties in this figure. For example, 
based on measurements of S 13 C in ocean sediments, Raymo et al. 
(1996) cite a mean value of 380 ppmv with maxima as high as 
425 ppmv. More recent data from Seki et al. (2010), using alkenones 
and boron isotope proxies, cite a mean of 360 ppmv with uncer- 
tainties + — 30 ppmv. Other recent data (Pagani et al., 2010) supports 
a mean of 380 ppmv. As such, for consistency with previous work, 
and to account for likely associated increases in non-C0 2 greenhouse 
gases such as are observed in the ice core record (Siegenthaler et al., 
2005), we consider here the effects of 350 and 450 ppmv as alterna- 
tive C0 2 concentrations. To first order, the temperature effects of ele- 
vated C0 2 are expected to scale logarithmically with the C0 2 
concentration. Therefore, it is possible to estimate the total mid- 
Pliocene temperature change for an arbitrary C0 2 level of x, AT C ° 2=X 
as: 


AT 


C0 2 =x 


d'l arog 3 ” dT veg 


+ ^ice + 


logjx/ 280) 
log( 400/280) c ° 2 


(15) 


For a C0 2 level of 350 ppm this gives AT C ° 2 = 350 = 2.7 °C, and a 
partitioning (see Table 1) of 1.0 °C (36%) from the C0 2 (dT C o 2 ), 
0.7 °C (26%) from the orography, 0.7 °C (26%) from the vegetation 
(dT V eg), and 0.3 °C (12%) from the ice sheets ( dT ice ). For a C0 2 level 
of 450 ppm this gives AT C ° 2 = 450 =3.8 °C, and a partitioning (see 
Table 1) of 2.1 °C (55%) from the C0 2 (dT C o 2 ), 0.7 °C (18%) from the 
orography, 0.7 °C (18%), from the vegetation (dT veg ), and 0.3 °C (9%) 
from the ice sheets ( dT ice ). 

Furthermore, given a ‘true’ mid-Pliocene global mean temperature 
change, A T c ° 2=x , we can solve Eq. (15) for x. By converting the 
PRISM3 estimates of global SST to estimates of global surface air tem- 
perature using a scaling factor, Lunt et al. (2010a) estimated the true 
AT C ° 2=X to be about 0.27 °C greater than the AT predicted by the 
model. This allows us to estimate x, the ‘true’ value of mid-Pliocene 
C0 2 , to be 380 ppmv. It should be noted that this calculation assumes 
that our uncertainty in C0 2 is much greater than uncertainties which 
arise due to model error, errors in the applied mid-Pliocene boundary 
conditions, and errors in the PR1SM3 SSTs. 



Fig. 6. Surface air temperature change due to C0 2 alone calculated as (a) Eq. (2) and (b) Eq. (3). The difference between the two approaches (equal to the synergy, S in Eq. 5) is 
shown in (c). 
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Fig. 7. (a) Difference in mid-Pliocene surface air temperature between our simulation and that of Haywood and Valdes (2004). (b) The same, but for the mid-Pliocene anomalies, 
mid-Pliocene minus pre-industrial. 


4.3. Climate variability through the mid-Pliocene 

The mid-Pliocene spans approximately 300,000 yrs, and, although 
relatively stable compared to the Quaternary, does display climate 
variability on orbital timescales (Lisiecki and Raymo, 2005), which 
can be interpreted as a series of glacials and interglacials (albeit 
much smaller in magnitude than those of the Quaternary). By 
combining high resolution mid-Pliocene oxygen isotope and Mg/Ca 
measurements, Dwyer and Chandler (2009) identified six sea level 
highstands during the mid-Pliocene of between 10 and 30 m above 
modern, and several lowstands, including Marine Isotope Stage KM2 
in the middle of the mid-Pliocene, estimated to be 40 m below mod- 
ern. However, we carry out a single simulation to represent this entire 
time period. 

For the orography, this is probably not an issue, as changes in 
orography occur over much longer timescales than orbital fluctua- 
tions. However, the orbit, C0 2 , ice sheets, and vegetation likely varied 
significantly through the mid-Pliocene. The orbital forcing in our sim- 
ulations is that of modern. At 65°N in June, the modern forcing is close 
to the average forcing of the mid-Pliocene, the difference being 
— 15 Win - 2 compared to a maximum difference of +50Wm” 2 
during the mid-Pliocene (Lunt et al., 2008). For C0 2 , we have used 
400 ppmv whereas the record of Raymo et al. (1996) varies between 
330 and 425 ppmv. For ice sheets, the PRISM2 reconstruction is 
characterised by a sea-level increase of 25 m compared to modern, 
whereas Dwyer and Chandler (2009) find variations in global sea- 
level of + — 25 m compared to modern, encompassing glacial/ 
interglacial variability. The PRISM3 SST evaluation dataset does con- 
sist of sub-orbitally dated sites. However, the PRISM SSTs do not 


Table 1 

Total mid-Pliocene global mean wanning compared to preindustrial (AT), and the 
global mean partitioning between C0 2 ( <ff C o 2 ), orography (dT orog ), vegetation ( dT veg ), 
and ice (dT,„). This is shown for the default case, and cases where the sensitivity to 
orography and C0 2 are tested, as described in Sections 4.1 and 4.2. 



AT 

ra 

dT C o 2 

ra 

dT or0 g 

[°C] 

dTi ce 

ra 

dT veg 

ra 

Default 

3.30 

1.58 

0.70 

0.70 

0.33 

Orography = modern 

2.60 

1.58 

0 

0.70 

0.33 

CO 2 = 350 ppmv 

2.71 

0.99 

0.70 

0.70 

0.33 

C0 2 = 450 ppmv 

3.83 

2.10 

0.70 

0.70 

0.33 


represent average SSTs through the mid-Pliocene but have been fil- 
tered via a process of ‘warm peak averaging' (Dowsett et al., 2009), 
which means that the PRISM3 SSTs represent average warm intergla- 
cial conditions in the mid-Pliocene. For vegetation, the data sites in 
the Thomson and Fleming reconstruction, upon which PRISM2 are 
based, are not dated to orbital timescale accuracy, and so each site 
could represent either glacial or interglacial-type conditions. The 
same is true of the Salzmann et al. (2008) vegetation dataset, with 
which our simulation has been evaluated. However, in locations 
where a number of possible biomisations were consistent with the 
data, Salzmann et al. (2008) chose the warmest, to maintain consis- 
tency with the SST warm peak averaging. 

As such, our simulations are a hybrid representation of the mid- 
Pliocene: the orbit, vegetation and orography being close to mid- 
Pliocene average, and the C0 2 and ice sheets being closer to intergla- 
cial values. The mid-Pliocene simulation has previously been 
compared with vegetation data which represent an average-to- 
warm mid-Pliocene palaeoenvironment (Lunt et al„ 2010a), and SST 
data which represent interglacial values (Dowsett et al., 2011). 
These discrepancies may go some way to explaining some of the 
model-data disagreements. For example, the greater high-latitude 
warmth in the PRISM SST reconstruction compared to the model 
could be a result of the warm-peak averaging, which by definition 
biases the SST reconstructions to warm values. Future work will aim 
to carry out simulations more representative of specific time periods 
within the mid-Pliocene, and to compare these to orbitally-resolved 
versions of the PRISM SST dataset. 


4.4. Model uncertainties 

Uncertainties associated with the model itself (as opposed to the 
boundary condition uncertainties discussed above) can be broadly 
divided into 'parametric uncertainty’ and 'structural uncertainty'. 

Parametric uncertainty relates to uncertainties in model parame- 
ters. These parameters are often associated with the representation 
of sub-gridscale processes and include, for example, the gridbox- 
average relative humidity at which clouds are assumed to start form- 
ing. They are generally poorly constrained by observations and so are 
essentially ‘tunable’. A single model simulation, as presented in this 
paper, can only represent one single point in the whole space of pos- 
sible plausible parameter combinations, and as such undersamples 
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the range of model possibilities. The full space can be explored by car- 
rying out simulations in which these tunable parameters are per- 
turbed. A preliminary study has been carried out with this model in 
the context of the mid-Pliocene (Pope et al., 2011). That study 
found a range of AT of 2.7 °C to 4.5 °C and could therefore be used 
to place approximate error bars on our AT; however, it did not inves- 
tigate the causes of such a change, so the impact of uncertain param- 
eters on our factorisation is unclear, and is a focus of ongoing work. 

Structural uncertainty relates to changes in the model which can 
not be made purely by modifying the values of tunable parameters. 
It relates to our uncertainty in the physical processes themselves 
which govern Earth System behaviour, and our inability to implement 
complex processes in a numerical model of a given resolution. Some 
information on the magnitude of this error can be obtained by consid- 
ering other climate models. Haywood et al. (2009) compared two 
structurally different models, of the mid-Pliocene. They had a range 
of AT of 2.39 °C to 2.41 °C (this is very much a minimum uncertainty 
range, especially as those simulations were carried out with 
atmosphere-only models). Again, whilst putting some context to 
our results, it is not clear how this uncertainty would affect our fac- 
torisation or energy balance analysis. Ongoing work, in the frame- 
work of the project PlioMIP, is aiming to gain more information on 
the structural uncertainty by comparing many atmosphere-only and 
atmosphere-ocean Pliocene simulations produced by different 
models (Haywood et al., 2010, 2011b). 

5. Conclusions 

Using a novel form of factorisation, we have partitioned the causes 
of mid-Pliocene warmth between C0 2 (36%-61%), orography 
(0-26%), vegetation (21%-27%) and ice sheets (9-13%). The ranges 
are estimated by considering the sensitivity of the results to uncer- 
tainties in the mid-Pliocene C0 2 concentration and orography (sum- 
marised in Table 1). Despite the relatively small contribution of ice 
sheets on a global scale, it is responsible for the majority of Southern 
Hemisphere high latitude warming. Northern Hemisphere high- 
latitude warming is due mainly to a combination of C0 2 and orogra- 
phy changes. Furthermore, we have carried out an energy balance 
analysis, and shown that surface albedo changes and direct 
greenhouse-gas forcing contribute significantly more than cloud 
feedbacks to the total mid-Pliocene warming, with the greenhouse 
gas forcing dominating in low latitudes, and the surface albedo 
changes dominating at mid-high latitudes. 

Future work should further assess the sensitivity of these results 
to the boundary conditions applied (for example by using the 
newer PRISM3 reconstructions compared with PR1SM2 used here, 
and extending the datasets to include varying soil properties), to 
the model used, and to parameters within the models themselves. 
Both the modelling and data communities should start to investigate 
orbital-scale variability within the mid-Pliocene. This is particularly 
important for assessing the real relevance of the mid-Pliocene as an 
analogue for long-term future (sub-orbital timescale) climate change. 
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